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We present a detailed derivation of a multi-site mean-field theory (MSMFT) used to describe 
the Mott-insulator to superfluid transition of bosonic atoms in optical lattices. The approach is 
based on partitioning the lattice into small clusters which are decoupled by means of a mean field 
^Sj . approximation. This approximation invokes local superfluid order parameters defined for each of the 

boundary sites of the cluster. The resulting MSMFT grand potential has a non-trivial topology as a 
c 2 ^ ■ function of the various order parameters. An understanding of this topology provides two different 

, criteria for the determination of the Mott insulator superfluid phase boundaries. We apply this 

formalism to d-dimensional hypercubic lattices in one, two and three dimensions, and demonstrate 
^ ■ the improvement in the estimation of the phase boundaries when MSMFT is utilized for increasingly 

larger clusters, with the best quantitative agreement found for d = 3. The MSMFT is then used 
fT^ . to examine a linear dimer chain in which the on-site energies within the dimer have an energy 

separation of A. This system has a complicated phase diagram within the parameter space of the 
model, with many distinct Mott phases separated by superfluid regions. 
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I. INTRODUCTION 

The properties of both the pure and disordered Bose-Hubbard (BH) model were first elucidated in a remarkably 



^4 insightful paper by Fisher et al. 1]. For the case of a pure (or homogeneous) system, it was shown that a collection 
• \ of interacting bosons on a lattice undergoes a Mott insulator (MI) to superfluid (SF) transition as the ratio J /U is 

■ varied, where J characterizes the rate of hopping from one lattice site to another, and U represents the interaction 
between two bosons residing on a given site. The insulating phase appears as lobes in the ^/U vs. J /U plane, where 
/i is the chemical potential, within which the average occupancy per site takes on an integral value. Outside these 

''O ' bounded regions, the system is a superfluid characterized by a nonzero condensate order parameter, and the density 
C . varies continuously with the system parameters. 

^ ■ The BH model became particularly germane with the advent of trapped Bose gases and Jaksch et al. [2] proposed 
I i ' that the model would provide a realistic theoretical description of bosons residing in an optical lattice. They obtained 
estimates of the parameters appearing in the BH model and argued that the model is relevant to the physical systems 
experimentally accessible. Subsequent experiments [1] indeed showed that trapped Bose gases are an ideal setting 
J^t . within which to study the theoretically predicted MI-SF transition. 

■ Many theoretical studies of the BH model followed the original Fisher et al. paper [l| using a variety of theoretical 
methods and approximations and the properties of the transition in one, two and three dimensions are now well 
established . Apart from the considerable work done on the disordered BH model [H 0) S d; [3 j the model 
has also been extended to su per lattices T, 14-22], spinor condensates [23 25], multi-component systems 32, 26 2^ 
and multiband situations |29l - [3l| . With increasing complexity, the MI-SF phase diagram becomes increasingly richer 
in structure. 

One of the most useful theoretical approaches for obtaining a qualitative understanding of the MI-SF transition 
^ ' is mean- field theory which was originall y rn otivated by considering the infinite-range hopping limit [l|. Subsequent 
, reformulations of mean- field theory [s^lssj invoked the existence of a condensate order parameter which was used 

■ to decouple the nonlocal hopping term of the BH Hamiltonian. In its simplest form, which we refer to as the site- 
H , decoupled mean- field theory (SDMFT), the decoupling leads to a system Hamiltonian consisting of a sum of site 

. . . ■ Hamiltonians. The latter are effectively independent but depend on the order parameter, in general site-dependent, 
which can be thought of as a variational parameter. For a homogeneous lattice, the ground state of the system is 
determined by minimizing the system energy with respect to the order parameter. If the energy is minimized for a 
non-zero value of the order parameter, the system is in the SF phase, otherwise it is in the MI phase. The phase 
boundaries obtained using this approach are consistent with the results of more sophisticated approaches. It can 
be shown [32| that SDMFT is equivalent to the alternative starting point based on the Gutzwiller ansatz for the 
ground-state wave function [H, Ugl . 

One of the limitations of the SDMFT is the neglect of inter-site correlations which allow for fluctuations of various 
physical variables. For example, in the MI phase of a homogeneous system, the site occupancy is precisely integral, 
whereas in reality some (albeit small) fluctuations in the site occupancy must occur. These effects can be captured, at 
least to some extent, by dividing the system into clusters of arbitrary size and using mean-field theory to decouple the 
clusters. We refer to theories of this kind as multi-site mean-field theories (MSMFT) and in this paper, explore this 
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approximation for various situations. This approach has been used previously by Buonsante et al. |15l . |17[ to study 
the MI-SF transition in superlattices where novel (loophole) features emerge and in an investigation of the disorderd 
BH model It should be noted that a two-site mean- field theory was also introduced in [s^l using the so-called 
phase-space method. 

The main purpose of this paper is to provide a detailed derivation of the MSMFT and to explain how it can be 
used to systematically improve the quality of the results for the MI-SF phase boundary. It should be emphasized at 
the outset that the method cannot compete in a quantitative sense with more sophisticated methods such as quantum 
Monte Carlo (QMC) [4, 7-9], the density- matrix renormalization-group (DMRG) method ^^.^ other theoretical 

techniques [a, [jLlI- However, the MSMFT method is relatively straightforward and allows one to explore efficiently 
the dependence of its predictions on the various parameters which define more complex physical models. In addition, 
as stated earlier, it has the merit of providing useful information about inter-site correlations that are missed in the 
SDMFT. 

Our paper is organized as follows. In Sec. Ill Al we derive the MSMFT for the simplest case of a one-dimensional 
lattice where the definition of clusters used in the multi-site decoupling is particularly straightforward. The method 
is extended in Sec. IIIBI to higher dimensions where more freedom is available, within certain limits, in defining the 
clusters which cover the lattice. The MI-SF transition is determined by means of a grand potential which is a function 
of the various superfluid order parameters defined for the cluster. A detailed analysis of the grand potential leads to 
two different criteria for the determination of the phase boundaries. Our results for the case of homogeneous lattices 
are presented in Sec. IIIII In Sec. IIVI we consider one-dimensional superlattices; the simplest consists of alternating A 
and B sites which we refer to as a dimer chain. In this case, the nonequivalence of the sites requires two different 
order parameters for the description of the superfiuid phase. We also consider an example of a four-site superlattice 
where qualitative differences appear in the predictions of SDMFT and MSMFT. 

II. FORMALISM: MULTI-SITE MEAN-FIELD THEORY 

A. Derivation in One Dimension 

Our work is based on the BH Hamiltonian '^l in the grand canonical ensemble. With the assumption of a single 
orbital per site, this Hamiltonian is given by 

iC^a - ^ = X! " ^^)ni + ^'^Uini{ni - l)-^jycjcj, (1) 

i i ij 

where the index i labels the sites of the optical lattice and cj and Ci are site creation and annihilation operators 
(henceforth referred to as site operators), respectively; the number operator for site i is given by fii = c[ci. The 
system parameters include the on-site energies at each lattice site, the tunnelling energy Jij between sites i and j, 
and the intra-site interaction energy Ui. To a good approximation it is sufficient to ignore interactions between bosons 
on different sites and hopping between sites further apart than the nearest-neighbour distance 0, [ssj . Furthermore, 
we will restrict our considerations to the case where the interaction parameter has a common value U for all sites 
and a hopping parameter J for all nearest-neighbour pairs. Generalizations to more complex situations such as 
superlattices jla - KLTj can be readily accommodated in the MSMFT that we develop. The final parameter in the BH 
Hamiltonian is the chemical potential /i which controls the number of particles in the system. Although extensions 
to finite temperatures are certainly feasible (l6l |. we will only consider the properties of the BH Hamiltonian at zero 
temperature. 

The MSMFT [l^ is most easily formulated for the example of a homogeneous one-dimensional chain with nearest- 
neighbour hopping. The on-site energies Si can be set to zero and the interactions are taken to be site independent. 
The first step in the derivation is to partition a chain having Ns sites into clusters each containing L sites, so that 

Ns = LN,. (2) 

We will refer to the cluster of length L as an "L-mer" . A schematic of the partitioning being considered is shown in 
Fig. [H The hopping part of the Hamiltonian, i-Lhop — ^JJ2{ij) ^l^j^ '^^^ then be written as 

i=0 1 = 1 3=0 
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FIG. 1: The partitioning of a one-dimensional chain into hnear clusters of length L. Within our mutli-site mean-field 
formulation each cluster with open boundary conditions is treated exactly, and the inter-cluster couplings (generated by the 
hopping Hamiltonian) are treated in a mean-field decoupling approximation. 



where we have isolated the terms with / = i in the second sum which couple sites between adjacent L-mers. This 
term will be denoted by Hcoup- We assume periodic boundary conditions so that cln +i = ci. 

The coupling between L-mers can be eliminated by invoking the usual argument [33 . [33j . We assume the existence 
of a homogenous superfluid order parameter 



(4) 



and write Ci = •ip+{ci — 'ip) for each of the operators in Hcoup- Neglecting quadratic terms in the fluctuation Sci = Ci — ip, 
we find 



J=0 



(5) 



We thus arrive at the cluster-decoupled Hamiltonian 



J=0 



(6) 



where JC^ only depends on the site operators within the j-th L-mer. Taking j — 0, we have 
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and 



^ = - J 



(7) 
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The V'- independent operator K.^ is the Hamiltonian of a one-dimensional chain of length L with open ends. The dif- 
ferent L-mers are independent physical systems but are effectively coupled by means of the order parameter appearing 
in the mean- field perturbation Vff^ . 

It is clear from the form of © that the phase of the order parameter, if] — 1-01 exp (it/)), can be absorbed by a 
redefinition of the site operators: exp {—i(j))ci — >■ c;. The resulting Hamiltonian depends on [ipl and we can therefore 
take the order parameter to be a real, positive quantity. The ground state of the system described by a single 
homogeneous order parameter can be determined by minimizing ICq^^ with respect to ip in th L-mer Fock space. 



B. MSMFT in Higher Dimensions 



The extension of the MSMFT to higher dimensions is straightforward, although some new elements appear due to 
the freedom available in partitioning the lattice into clusters. This is illustrated in Fig. 2(a) for the example of a 
two-dimensional square lattice. The figure shows how the lattice can be covered using L-mers with size L = 1, 2, 4, 5 
and 6. It is clearly necessary that these clusters cover the entire lattice without duplication. It is also desirable that 
they have the same point-group symmetry as the original lattice. The examples L = 1, 4 and 9 satisfy this latter 
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FIG. 2: The partitioning of a two-dimensional square lattice using various repeated clusters. Shown in (a) are single site 
monomers, dimers (L = 2), and L-mers for L = 4, 5, 6. Shown in (b) is a portion of a two-dimensional square lattice that is 
partitioned into 3x3 clusters. The solid lines connect sites in a given cluster, dotted lines show the division of the lattice into 
such clusters, and the dashed lines show the inter-cluster couplings that are treated within mean-field theory. In each cluster 
there are three inequivalent sites: the central site C, which does not couple to neighbouring clusters, and corner and edge sites, 
which respectively are labelled A and B. 



criterion. The L = 5 cluster has the point-group symmetry of the square lattice, but its covering of the 2D plane 
introduces a chirality not present in the homogeneous system. 

Although we develop the MSMFT in its most general form, it will be useful to keep in mind the example of a 
two-dimensional square lattice with the 3x3 clusters shown in Fig. 2(b). The sites in the lattice can be specified by 
R-f T, where R is a cluster Br avals lattice vector and r defines the basis of sites within the cluster. For the example 
of Fig. 2(b), r runs over the nine sites of the 3x3 clusters. A site on the boundary of a cluster will be denoted by a 
and is connected, via J, to one or more boundary sites on adjacent clusters. The mean-field decoupling of the hopping 
term for the boundary site a in the L-mer and the boundary site /?' of the adjoining L'-mer is achieved using the 
prescription 

In general, we allow for different order parameters for each of the boundary sites. Since the (3' site in th L'-mer is 
related to the /? site in the L-mer by translational symmetry, we have denoted the order parameter on this site as 
Tpp. (For example, the sites A'l and Ai in Fig. 2(b) must have the same order parameter by translational symmetry.) 
For now, we have also allowed the order parameters to be complex. Collecting all terms pertaining to the L-mer, the 
mean-field decoupling leads to the L-mer mean-field perturbation 

vf ^({^ J) ^-J2jc.p {ci^p + c^rp - cv'^) , (11) 

where the sums extend over all boundary sites. To obtain this form we have defined the symmetric matrix j7 with 

matrix elements Jap = J gap where gap is the number of clusters to which the L-mer of interest is connected by a 
pair of a and /? sites in the way described above. If such a pair of sites is not coupled, g^p = 0. We will refer to J_ as 

the connectivity matrix, since it enca psu lates the way in which the different sites in the cluster are connected to one 
another via the mean-field couplings [34|. 

The perturbation in (|lip depends on the set of order parameters {il^a} defined for the ensemble of boundary sites. 
The cluster Hamiltonian then takes the general form 



(12) 
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where the ■0Q,-independent part 
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(13) 



is the suni of terms in the grand canonical Hamiltonian K, that depend only on the L-mer variables. The final sum in 
(1131) is restricted to nearest-neighbour pairs within the L-mer. It is important to note that this intra-cluster hopping 
leads to inter-site correlations even within the Mott insulating phase. 

The MSMFT Hamiltonian depends on various factors, such as the dimension and geometry of the lattice, the shape 
and size of the clusters and the coordination of sites in one cluster with those of its neighbours. The extent to which 
different order parameters are required depends on the physical application and will be clarified by example. However, 
we emphasize that these order parameters are not prescribed but in general are determined by the solution of the 
mean-field problem itself. 

The cluster grand canonical Hamiltonian (jl2p is Hermitian but it is not number-conserving due to the perturbation 
in (fTTj) . As a result, its eigenvectors must be determined in the i-mer Fock space. We denote the state with the 
lowest eigenvalue ^oiii^a}) as |5'o)- For a normalized state vector, the variation of r2o({V'a}) with respect to ■(/'^ is 
given by 



*0 



9^; 



This implies that rio({^Q}) is stationary when the order parameters satisfy 

{C^) = (*o({V^a})|c^|*o({^a})) = 4^-, 



(14) 



(15) 



for all boundary sites 7, where the values of the order parameters at the stationary point are denoted by ^"7- We 
thus see that stationarity of the grand potential is associated with the physically necessary condition that the order 
parameters are determined self-consistently. If more than one stationary point arises, the physical state of the system 
is assumed to correspond to the stationary point with the minimum value of f2o({V'a})- 

Intuition might lead one to expect that a stationary point is an extremum of the grand potential, but it is straight- 
forward to show that it is not. Writing t/jq — i-'a + a and expanding (jl2l) about a stationary point, we have 



^''^({V^a}) = /C*^^({^J) - [(cl - C)AV'/3 + (Ca - Via) " ^Vc^Mp] 



(16) 



We denote the sum by I^Vff^ and consider it as a perturbation to the grand Hamiltonian /C^^^({?/'q,}) at the stationary 
point. The first order correction to the grand potential ^o{{i^a}) is 

1]W({AV;J) 



(17) 



while the second order correction is 



'^iEcbJcp ciA^p + c^Ar^ |o)i 



(18) 



where the states {v) are eigenstates of X^^^-'' {{ipa}) with eigenvalues n^dipa})- ^'o^ is negative definite and fl^^' can 
always be made negative by choosing Aipp = —Aipa for some pair of deviations with all others equal to zero. Thus, 
unlike the situation for the case of a single order parameter, a stationary point is not a local minimum; the value of 
^o{{'4'a}) can always be made smaller than r2o({V'a}) by moving away from the stationary point in some direction. 
For multiple order parameters, the stationary point is, in general, a saddle point and as a result, it cannot be located 
by means of a variational principle. Below, we provide criteria for identifying the emergence of a superfluid phase 
from a Mott insulating phase without having to appeal to a variational principle. 



C. Perturbative Treatment of the MI-SF Transition 



The point {ipa} = {0} is always a stationary point and it too is not an extremum in general. This point has special 
significance in that it corresponds to the Mott insulating phase. Assuming the MI-SF transition to be continuous 
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as a function of the system parameters (for example, the hopping strength J), we expect the stationary point to 
move continuously away from the {ipa} = {0} point as the parameters are varied beyond some critical values. This 
behaviour can be analyzed by treating (fTTI) as a perturbation to ■ Depending on the values of the system parameters 
(/X, J and U), /C° has a ground state |0) containing N particles. The range of parameters for which this is the case 
defines what we refer to as A^-domains in the multi-dimensional parameter space. Each of these A^-domains can, in 
principle, give rise to a Mott phase with a certain number of particles per cluster. These regions are conventionally 
referred to as Mott lobes. 

The first and second order corrections to the energy r2o({0}) are obtained with the replacement 'ipa ^ and 
Aipa i^a in (HZl and (HH]). We have 

ni'H{4,a}) = J2J^p^*^'Pp (19) 

and 

^ I ('^i + c^rf\ |0) p 

where in this case, the states \v) are eigenstates of /C° and are therefore number eigenstates. For the matrix elements 
in the sum to be nonzero, the state \v) must have either — 1 or A'^ + 1 particles. Defining the operator 

Oq = ^ JaP C/3, (21) 

we see that the second order correction is given by 

a/3 

where the Hermitian matrix M is defined by 

{o\daW){iy\dl\o) + {o\dlw){iy\6a\o) 

In obtaining this result we have used the fact that the operator Oa only has finite matrix elements between states 
whose particle numbers differ by one. In view of (pTj) . we see that 

M = J (24) 
where the elements of the Hermitian matrix C are given by 

{0\caW}{,.\cl\0} + {0\cl\,.}{,.\ca\0) 

"'"^0 mo})-n4{o}) • ^''^ 

Combining (ITOl) and (f^ . we have 

fio({^«}) = ^^o({0}) + ^ Wap CV'^ + • • • , (26) 

a/3 

where 

W = 1 + M ^ 1 + IQl. (27) 

We will refer to the Hermitian matrix W as the energy matrix. 

If the value f2o({0}) is the minimum of all stationary points, one is in the Mott insulating phase. The question 
then arises as to whether the stability of the Mott phase can be established independently of determining f2o({''/'a}) 
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for all of these points. Since J7o({0}) is not an extremum, it is not obvious a priori what properties the energy matrix 
W must have in order for the Mott phase to be stable. To address this query we determine the eigenvectors and 

eigenvalues of W: 

W■V^=iO,■V^. (28) 

Since W is Hermitian, the eigenvalues are real and the eigenvectors can be chosen to be an orthonormal set. Expanding 
the order-parameter vector tp = {t/ja} as 

V'^E^^^^" (29) 

i 

we see that 

Ar!oW=E^^I^^I'- (30) 

i 

The eigenvectors thus define directions in the order-parameter space along which the grand potential varies quadrat- 
ically with a curvature that is determined by the sign of uji. Furthermore, as one moves along a line in one of these 
directions, the grand potential is stationary with respect to displacements away from the line. Thus the transition to 
the SF phase must be accompanied by the appearance of a stationary point with xp ^ along one of these directions. 
The only scenario consistent with the continuous movement of the stationary point away from ■0 = is that one 
of the positive eigenvalues passes through zero. The first positive eigenvalue to do so establishes the criterion for 
the transition to the SF phase. This is very analogous to the situation one finds with a single order parameter [ssj . 
We elaborate on these properties of the grand potential below, and illustrate the general behaviour with a concrete 
numerical example in Sec. Ill El 

We now denote the vector relevant to the transition as Vi. When wi passes through zero, the stationary point moves 
continuously from xp = to a point where the order parameter is non-zero. If higher order terms in the perturbation 
expansion are retained, one expects the grand potential to behave, in the vicinity of the transition, according to the 
Landau theory [39} expression 

AfioW = ^|Ci|' + |lCi|', (31) 
where r — 2uji, and we assume u to be positive. For small and negative Wi, this leads to the stationary order parameter 

(32) 

Since the vector vi can be chosen to be a real, one sees that, in the vicinity of the phase boundary, the order 
parameter is real. In other words, there are no phase differences (apart from a possible sign) between the order 
parameter components. In fact, in all the examples we have studied, the components of the order parameter vector 
all have the same sign, which we take to be positive. One might argue that this is the expected behaviour, since 
relative sign differences in the components would lead to more rapid spatial variations of the order parameter with a 
resultant higher kinetic energy. However, as stated earlier, the grand potential is not an extremum at the stationary 
point and can take on lower values if one moves away from the stationary point in some direction v.; . Since this vector 
is orthogonal to vi, it must necessarily have components with different signs. One therefore cannot argue that the 
components of the physically relevant direction (vi) must have the same sign on the basis of energy considerations. 
Henceforth, we assume that the components of the order parameter are real. This assumption is supported by all 
self-consistent solutions of (|12l) that we have obtained in the SF phase. 

Referring to Eq. ((27|) , we see that the eigenvalue spectrum of W is closely related to that of the connectivity matrix 

J. In fact, for J -> 0, the eigenvectors and eigenvalues of these two matrices coincide. The eigenvectors of are 

defined by 

Juj = XiUi. (33) 

The eigenvectors with zero eigenvalue play a special role in that they are simultaneously zero-eigenvalue eigenvectors 
of W. These vectors are independent of the magnitude of J and other system parameters, and are simply determined 
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by the structure of J. The eigenvectors of W with non-zero eigenvalues can be determined by using as the basis 
vectors. Then, in what we cah the J-representation, we have 



(34) 



If we order the Xi = eigenvectors as i = 1, . . . , m and the remaining n eigenvectors as i = m + 1, . . . ,m + n, the 
non-zero eigenvectors of W are found in the rt-dimcnsional subspace spanned by the eigenvectors of with non-zero 

eigenvalues, and are determined by 



Here we have defined the n x n matrix 



where Dnxn is the diagonal matrix 



det (wyj„ ~ ^/„> 

■ \ 



(35) 
(36) 



( Xm+\ 

V ■■■ 





A„j+„ 1 



(37) 



and 



(J) 



m+l,m-j-l 



(J) 

m+l,m+n 



7n-{-n,m-\-n 



(38) 



The MI-SF transition is located by following the eigenvalues determined by ([HS)) as a function of a system parameter, 
such as J. As discussed above, the positive eigenvalue uji which first passes through zero determines the phase 
boundary. In view of ([55|) . the condition for this to happen is det(W^xn) = 0. 

Although the eigenvalue problem in psp eliminates the eigenvectors with zero eigenvalues, in practice it is more 
straightforward to determine all the eigenvalues and eigenvectors using When the uji eigenvalue is identified 

and found to pass through zero, the corresponding eigenvector directly determines the relative magnitude of the order 
parameter components just as one enters the SF phase. The main advantage of the J-representation is that is reveals 
the mathematical structure of the W matrix and facilitates some of the formal developments that follow. 

An alternate approach to that of following the W eigenvalues is based on constructing a stability criterion for the 

Mott phase [13] ■ It too is based on a perturbative analysis and can be derived as follows. To first order in perturbation 
theory, one finds that 



(39) 



where the matrix which we refer to as the stability matrix, is defined in terms of matrices introduced earlier, 
namely 

s = -c; J. (40) 

We note that this matrix is not Hermitian. Since the left-hand side of Eq. ([5^ is the perturbative estimate of ipa, 
this equation is suggestive of an iterative scheme defined by the linear map 



(fe+1) 
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where k = 0, 1, 2, . . . can be thought of as an iteration index with k 
Defining the eigenvectors of 5 by 



(41) 



oo identifying the self-consistent solution. 



(42) 
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we see that the iterative sequence converges to the Mott insulating = fixed point if the absolute values of the 
eigenvalues tJ^ are all less than one. Therefore, this procedure provides an operational definition of the Mott phase in 
terms of the eigenvalues of the stability matrix jS. 

The definition of the stability matrix in ()40|) shows that the eigenvectors of with zero eigenvalues are also zero- 
eigenvalue eigenvectors of S^, that is, z,j = for i = 1, . . . ,m, with ai = 0. In the J-representation, the stability 
matrix thus has the block structure 

giJ) ^ f Omxrn S^^l^ \ ^ ^^^^ 
\ Onxm "Snxn / 

where the null matrices (O) have the dimensions indicated, 

^nxn = ~^nxn^rixn (44) 

and 



X n 



(45) 



It is clear from that the cr = eigenvalue is m-fold degenerate and that the remaining eigenvalues of are 
determined by the equation 



det 5, 



(5ii)„ - a/„xn) = 0. (46) 



If these latter eigenvalues all have a magnitude less than one, then i/? = is a fixed point of the map defined in (|4ip , 
and the Mott phase is stable. 

The complete set of eigenvectors of consist of the m zero-eigenvalue eigenvectors of j7 together with the remaining 

n eigenvectors determined by considering the n-dimensional eigenvalue problem 

s'^Ly.^o^Y.. (47) 
In view of (|44l) . this equation is equivalent to the generalized eigenvalue problem 

CnLy\ = -0•^^«xny^, (48) 

where y- = Dnxnji- Since C^x^„ is Hermitian and £>~xn (the inverse of ([571) ) is a real diagonal matrix, it is easy 
to show that the eigenvalues ai are real. In addition, for distinct eigenvalues, we have the orthogonality relation 
{Dnxn yiY'yj — 0. Once the vectors have been determined, the S}'''^ eigenvectors for i = m + 1, . . . , m + n are 

given by 



(49) 



with X, = cr^ ^5'^Ly»- 

We now establish the connection between the two different criteria derived above. The energy and stability matrices 
are related by 

lil - S). (50) 

An explicit expression for S_ cannot be obtained from this equation since ^7 does not in general have an inverse. 
However, within the J-representation, we have the more useful relation 

^nxn ~ Dnxn{Inxn — S^^xn)- (^1) 
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Since Dnxn has an inverse, we can write 



'-'nxn ^nxn ^ nxn " nxn- 



D'^ W 

^-^ ^ V -») y ' r, 



iJ) 



(52) 



From this we see that an eigenvector 1, • • • ,qm+n)'^ of W^'^li ^i^h zero eigenvalue is an eigenvector of s'^^^^ 
with eigenvalue one. This particular zero eigenvalue of the W matrix is obtained at the critical values of the system 

parameters and the corresponding eigenvector, previously denoted by vi, is given by Vi — li'^i- Using the 

structure of the S}"''' eigenvectors in we thus see that 

S;vi = vi, when Wv^ = 0. (53) 

We have thus proved that the energy and stability criteria for the transition to the SF phase are equivalent and can 
therefore be used interchangeably. 

D. Self-Consistent Superfluid Solutions 

The previous subsection outlined two different algorithms for identifying the parameter values of the Bose-Hubbard 
model for which the MI phase is stable, and an example of the application of these methods is presented in the next 
subsection. When the MI phase loses stability, the T = ground state of the Bose-Hubbard model corresponds to a 
superfluid state. The order parameter is then determined by finding the ground state solution of (|12p which satisfies 
the self-consistency condition given by (ITSl) . 

It is a simple matter to obtain such self-consistent solutions numerically. The iterative map 



V-r^) - (*o(Wi'='})|fi.l*o({^i'=^})) 



(54) 



provides updated order parameter components in terms of their current values. We find that any initial guess of the 
order parameters, provided that they all have the same sign, converges to a unique self-consistent solution without 
the need for special mixing [40| or more advanced techniques. This iterative method was used extensively in a recent 
paper [isj on disordered systems where a large number of sites with different order parameters are required. For 
homogeneous systems with only one order parameter, Brent's algorithm for root finding is particularly successful ^l|. 

E. Example Calculation 

We now illustrate the calculation of the phase boundary and superfluid order parameters using the methods outlined 
in the previous two subsections. It is instructive to consider the 2D square lattice partitioned into 3x2 clusters with 
the site labelling indicated in Fig. [31 For this L-mer, there are six boundary sites: the four corner (A) sites and the 
two edge (B) sites. The choice of this particular cluster is useful for a number of reasons. First, it demonstrates 
that meaningful physical results can be obtained even though the cluster does not have the symmetry of the two- 
dimensional square lattice. Second, it is a relatively simple example with several order parameters which illustrates 
the general formalism and the structure of the solutions. And lastly, it has obvious symmetries which allow us to 
reduce the number of independent order parameters; here there end up being two, one for the corner sites and one 
for the edge sites. Thus, for demonstration purposes this is an ideal cluster to examine. 



A4 




A3 












A] 


Bl 


Ai 





FIG. 3: The partitioning of tlie two-dimensional square lattice into 3x2 clusters. The cluster consists of inequivalent corner 
{A) and edge [B) sites. 
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The calculation of the MI-SF phase boundary can be approached in terms of the energy or stability criterion, each 
of which involves the connectivity matrix J^. For the site labelling shown in Fig. [3] and the order parameter specified 

as V' = (V'Ai , "0^2 7 , ''Pa4 , ^'Bi 1 4'B2 y I the connectivity matrix is given by 





J 





J 







J 





J 














J 





J 








J 





J 


























J 


V 
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0/ 



(55) 



We observe that the matrix is block diagonal reflecting the fact that corner sites are coupled to corner sites, and 
edge sites are coupled to edge sites. Following the notation of (pS)) . the j/ matrix has two A = eigenvalues, and the 

non-zero eigenvalues are ± J and ±2 J. The energy ^7} and stability pp]) matrices also depend on the C matrix in 

(p5| . The evaluation of this matrix requires the eigenstates and corresponding eigenvalues r2i/({0}) of the cluster 
Hamiltonian in (|13p with on-site energies ei — 0. We present details of the evaluation of these states in an Appendix. 

As a concrete numerical example, we choose a chemical potential of fj,/U = 0.4 which, for small J/U, places the 
system in the first Mott lobe with (n) — 1. The eigenvalues of the energy matrix are obtained by solving (|28p . 
whereas the eigenvalues of the stability matrix are obtained from (|42p . The four non-zero eigenvalues of these two 
matrices are shown in Fig. 2] as a function of J/U. We observe that one eigenvalue of the energy matrix associated 
with the eigenvector vi increases from zero to positive values, and then passes through zero at the critical hopping 
Jc/U = 0.04815 which locates the MI-SF phase boundary for this value of fJ,/U. At this same hopping one finds that one 
of the eigenvalues of the stability matrix first attains a magnitude of one, consistent with the discussion following (j52p . 
Figure [5] displays the variation of the eigenvector components of these two specific eigenstates, namely, the eigenvector 
Vi of (P5|) and the eigenvector Zi of ((i^ . Both of these eigenvectors have the form (1/1^, V'a, V'A, V'A, "05, '/'s)"^- One 
sees that these two eigenvectors are the same (zi = vi) at the MI-SF phase boundary and that (|53|) is indeed satisfied. 

Although we allowed for a six-component order parameter in the above analysis, we find that the relevant eigenvec- 
tors Vi and Zi only have two independent components corresponding to the two inequivalent sites in the 3x2 cluster. 
One, ^A, is associated with the four corner sites and the other, ■0b, is associated with the two edge sites. As can 
be seen in Fig. [5l these two values are almost equal for J/U close to the MI-SF transition. Even though the cluster 
does not mirror the symmetry of the square lattice, the violation of homogeneity is relatively weak. This behaviour 
persists into the SF phase. Fig. IH] shows that the self-consistently determined order parameters are close to each other 
over the range of hoppings indicated. The variation shown in Fig. [5] is consistent with the mean-field prediction in 
((32)) . as is the ratio of the two components. Also shown in this figure is the grand potential r2o({'0a}) for the order 
parameter = {tpA,'4'A,'4'A,'4'A,4'B,4'B)'^ as a function of ipA, where the ratio iIjb/4'a is fixed and taken from Fig.[S] 
at J/U = 0.04815. In the SF phase the grand potential shows minima at points which are very close to the order 
parameters determined self-consistently according to (jlSp . Thus, near the phase boundary, the grand potential is 
indeed given to a good approximation by the Landau expansion in ((3T|) . 

In view of the fact that the relevant eigenvector Vi only has two independent order parameter components, namely 
TpA and "05, it is of interest to see to what extent the calculations can be simplified by assuming an order parameter 
having the form tp = {ipA,ipA,tpA,ipA,'4'B,'4'B)'^ ■ With this assumed form, the eigenvalue problem in reduces to 

where Na and Nb are the number of occurrences of ^a and ^b, respectively, in the order parameter vector. For this 
specific example, Na = 4 and Nb = 2. The reduced matrix elements in (j56p are defined as 

^AA^J2J2^^i^ wab^J2J2^^i^ wba^J2J2^^^^ wbb^Y.J2^^^- (^7) 

i€Aj£A i£Aj£B ieB j£A i£B j£B 

By the same token, the expansion of the grand potential in (j26p is given by 

noi^AM = Oo(o,o) + ( 0. 0B ) ( Ztl){t)- ^''^ 

We note that the eigenvalues of the matrix W^°'^ appearing in this expression do not in general yield the desired W 
eigenvalues which must be obtained from the generalized eigenvalue problem in (j56l) . 
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FIG. 4: The non-zero eigenvalues of the energy (top 
panel, in units of U) and stability (bottom panel) ma- 
trices for a homogeneous two-dimensional square lattice 
using a 3 X 2 cluster vs. J/U. The chemical poten- 
tial is n/U = 0.4. One eigenvalue of the energy ma- 
trix crosses zero at J/U = 0.04815, the same position 
at which one eigenvalue of the stability matrix first at- 
tains a magnitude of one (the solid black curve in each 
figure). The superfluid phase is the T = ground state 
for J/U > 0.04815. 



FIG. 5: The normalized eigenvector components cor- 
responding to the energy (top) and stability (bottom) 
eigenvalues shown by the the solid black curves in Fig. |4] 
vs. J/U. The bracketed numbers give the number of 
components having the indicated values, and correspond 
to the corner sites (multiplicity 4) and edge sites (multi- 
plicity 2) of the 3x2 cluster. Note that at the transition 
these eigenvectors are the same, as required by ((53)) . 



Performing the sums in ((57)) . we obtain the explicit expressions 

w,, = Z.N.J + J g nom)-^.m) ' ^''^ 

WsB = zbNsJ + ZsJ E »,({0})-^,({0}) ' (60) 

= 5 ^om)~-^.m) ■ ^''^ 

Here we have introduced the notation (6^) for the site operators of the A (B) sites and define the connectivity 
number z. (zb)- For the specific example of the 3x2 cluster, we have = '2 and zb = 1- The element Wba is 



and 
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FIG. 6: The upper figure shows the self-consistently determined order parameters corresponding to the corner (iI'a) and edge 
(■0s) sites of the 3x2 cluster as a function of J/U. The arrows indicate the values of ipA at J/U = 0.049 and 0.05. The 
lower figure shows the grand potential as a function of tpA for the fixed ratio ipB/'tpA = 0.9396 found at the phase boundary at 
J/U = 0.04815. The curves from top to bottom are for values of J/U from 0.046 to 0.050 in steps of 0.001. In the SF phase, the 
grand potential shows minima at positions which are very close to the values of 'ipA found from the self-consistent calculations 
(indicated by the arrows). These results were obtained for n/U — 0.4. 

obtained from Wab by interchanging A and B. Since the matrix elements can be chosen to be real, we in fact have 
Wba = Was- 

The assumption of a two-component order parameter can also be used at the very beginning of this analysis to 
simplify the mean- field perturbation in ()11|) . One finds for the example being considered that 

V3*x2 = -ZAJi'Aj^ial + a,;) - zbJ^Jb Y^Cbl + hi) + NazaJ^^I + NbzbJ^^I- (62) 

A perturbative treatment of this mean-field perturbation will of course yield precisely the grand potential in (|58l) . 
However, the direct calculation of ([58]) does not by itself reveal the correct eigenvalue equation (|56|) required in the 
determination of the W eigenvalues. 

The eigenvalues of the energy matrix in (j56p are 



^± = liWAA + Wbb) ± \j ' + wl^^ (63) 
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where Wpq = Wpq / y/NpNg . In the Umit J — 0, we find that w+ ~ zaJ and aj_ ~ ZbJ- Thus, when symmetry is 
used to reduce the order parameter to the two components ipA and ipB,^Q see that the two energy matrix eigenvalues 
provided by the reduction increase hnearly with J and are positive. In fact, these two eigenvalues correspond to the 
two eigenvalues in Fig. IHlwhich increase positively from zero. We thus see that the reduced grand potential has a local 
minimum for small J. This behaviour should be contrasted with the result obtained using the full six-component 
order parameter for which negative eigenvalues of the energy matrix are found even in the limit of small J (see Fig. [6]). 

To determine the MI-SF phase boundary with increasing J we must find a zero crossing of one of the u)± eigenvalues. 
Since w+ > a;_ for all values of J, the eigenvalue which first goes to zero is w_ . Thus the location of the phase boundary 
is given by the condition 



which is equivalent to 



\{Waa + Wbb) - W \wls= 0, (64) 



WaaWbb = WIb. (65) 



We will use this equation in the following section to map out the MI-SF phase boundary for the two-dimensional 
square lattice. We note that the Na and Nb factors have dropped out in (|65|) as they must since they play no role in 
([55)1 when cj = 0. In other words, it is sufficient to look for the zero eigenvalue of W^^'^ in the determination of the 

phase boundary. Also, the ratio of the order parameter components just as one enters the SF region will be given by 



^A V WaA ■ ^ ' 

In light of the very similar numerical values of ipA and ipB shown above, it is reasonable to try one final simplification, 
namely, enforcing the identity ipA = ipB- When such an order parameter vector is used, the energy and stability 
matrices reduce to scalars and one recovers a single order-parameter theory. The critical hopping parameter obtained 
in this case for fi/U = 0.4 is Jc/U = 0.04816, as compared to Jc/U = 0.04815 obtained by allowing ipA and ipB to be 
different. Clearly, the inequivalence of different boundary sites in the MSMFT does not in itself adversely affect the 
prediction of the location of the MI-SF phase boundary for a homogeneous lattice. 



III. APPLICATION TO HOMOGENEOUS LATTICES IN 1, 2, AND 3 DIMENSIONS 

In this section we present our results for bosons on homogeneous lattices in d dimensions, focussing on the linear 
chain, square, and simple cubic lattices. We discuss the phase diagrams that are predicted by MSMFT, and the 
spatial correlations that are incorporated into ground-state wave functions found using this approach. 

First, we consider the phase diagram of the linear chain in the region of the first Mott lobe (0 < li/U < 1). 
The relevant MSMFT data is shown in Fig. [71 "Exact" numerical results, found from a DMRG study [6|, are also 
shown for comparison. The SDMFT result is seen to be far from the exact phase boundary, a result well known in the 
literature 6, 32] . In Fig. [7] we also show the results from the application of MSMFT for L-mers of size L = 2, 4, 8. One 
can see that as L increases the extent of the Mott lobe approaches the DMRG results. Therefore, in one dimension we 
find that there is a systematic improvement in identifying the phase boundary upon going from SDMFT to MSMFT. 
However, the predicted phase boundaries are still not close to the exact results and lack the cusp-like feature at the tip 
of the lobe. Further improvements could be achieved by increasing L but the progression in Fig. [7] suggests that going 
to L = 16 would provide only a modest improvement. This value of L is already beyond our numerical capabilities. 

Also shown in Fig. [7] are results obtained from an exact diagonalization of the full Hamiltonian in ([1]) for a finite 
linear chain with periodic boundary conditions. Since the Hamiltonian commutes with the total number operator, 
one can determine the iV-domains in the ji/U-J /U plane where the ground state of /C has a total of N particles. For 
chains of length A^g, the first Mott lobe must be found in the region where = A^^, at least in the limit Ns ^ oo. 
The boundaries of the N = Ns region are shown in Fig. [7] for chains of length Ng = 5 and 10. It is clear from the 
figure that the N = Ns domain bounds the exact Mott lobe with increasing accuracy with increasing A^^. However, 
exact diagonalizations for much longer chains would be necessary to accurately represent the tip of the Mott lobe. 

Our results in two dimensions are shown in Fig. [H In the MSMFT calculations we have used clusters of sizes 2x2 
and 3x3; the "exact" results are taken from the Monte Carlo study of Ref. As in the previous figure we also 
show the boundaries of the A^-domain (A^ = 9 in this case) for the exact solution of the full Hamiltonian of a 3 x 3 
cluster with periodic boundary conditions. To perform the calculations for the largest cluster it was necessary to take 
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FIG. 7: [Colour online] The MI-SF phase boundaries for a linear chain predicted by MSMFT, compared to the DMRG results 
of shown as circles. The phase boundary systematically approaches the exact result with increasing length L of the cluster: 
L — 2 (dash-dash-dot red line), L = 4 (dot-dash green line), L = 8 (dot-dot-dash blue line). Also shown are the boundaries of 
the A^-domain where the ground state of the full Hamiltonian for a finite chain of length Ns has N = Ns particles. 

advantage of the symmetry of the cluster; some details of the method used are given in the Appendix. Similar results 
are shown in Fig. [S] for three dimensions, where we have performed calculations for 2x2x1 and a 2 x 2 x 2 clusters. 
One sees a progressive improvement of the phase boundary with increasing cluster size, as in the case of the linear 
chain, with the phase boundary approaching the exact (MC) results The best agreement is obtained in three 
dimensions as one might expect on the basis of the validity of mean field theory when the upper critical dimension is 
approached. 

For the 3x3 cluster in Fig. HJb), the A and B sites are inequivalent and as a result, the order parameters ipA 
and ipB as one enters the SF phase will be different. We find, however, that the ratio given in (l66l) deviates from 
unity by no more than a few percent. It is therefore of interest to see the effect of enforcing the equality of ipA and 
tpB in the MSMFT calculations. (Note that this does not reduce MSMFT to SDMFT since one must still determine 
the eigenstates of /C" for the whole cluster.) We find that the phase boundaries shown in the two-dimensional phase 
diagram are virtually unchanged when this constraint on the order parameters is imposed. Thus the calculations for 
homogeneous lattices can be simplified by assuming that the order parameters of all the boundary sites of the cluster 
are the same. 

With several sites in the cluster, one can also consider spatial correlations between the various sites. As mentioned 
in the introduction, this is one of the main advantages of using a MSMFT. The correlations we consider are those 
exhibited by the function 

= {4>Uj), (67) 

where 

h = ai~ (a^). (68) 

In the MI phase, C'ij reduces to the single-particle density matrix while in the SF phase it is the contribution to the 
density matrix from the noncondcnsed component. 

We first consider our results for the linear chain, specifically along the fi/U — 0.4 line in the phase diagram of 
Fig. [71 For a given L-mer, we calculate Ci,i+i for all possible values of i and I consistent with the size of the L-mer. 
In the MI phase we find that Ci^i+i has a rather weak dependence on the position i within the cluster despite the fact 
that the states of the cluster are calculated with open-chain boundary conditions. For this reason we plot only Ci{L), 
the average of Ci^i+i over i. This quantity is shown in Fig. 1101 for L = 7 for I between and 6 and for values of J/U 
that span the phase boundary. Since this log-linear plot reveals a dependence on I which deviates weakly from a pure 
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FIG. 8; [Colour online] The MI-SF phase boundaries for a two-dimensional square lattice, as calculated using SDMFT (solid 
black line) and MSMFT (2 x 2, dashed red line; 3x3 (dotted blue line). The MC data (solid black circles) are taken from 
Ref. [^. Also shown (dot-dash black line) is the A'^-domain (A'^ = 9) for the 3x3 cluster determined by calculating the ground 
state of the full (non-MFT) Hamiltonian with periodic boundary conditions. 




FIG. 9; [Colour online] The MI-SF phase boundaries for a three-dimensional cubic lattice, as calculated using SDMFT (solid 
black line) and MSMFT (2x2x1, dashed red line; 2x2x2, dotted blue line). The MC data (solid black circles) taken from 
Ref. 



exponential, we have fit this data to the form 



where the fit parameters Iq, ^ and r] are in general functions of L [42]. We emphasize that this form of fitting function 
is chosen simply because it describes the short-range behaviour of the correlation function with reasonable accuracy. 
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FIG. 10: The average correlation function of the linear chain within MSMFT versus the intersite spacing. Panel (a) shows 
results for L = 7 in the MI phase: the six curves starting from 1 at n = correspond, from bottom to top, to J/U = 0.02, 0.04, 
0.06, 0.08, 0.10 and 0.12. The points are the calculated values of Ci{L) and the dashed curves are plots of the fitting function 
(|69|l . The lowest curve in (a) shows results for J/U = 0.10 for L = 3, 4, 5, 6 and 7, all plotted on top of each other. This curve 
has been displaced from the others for clarity by multiplying the data by 10~*. Panel (b) shows results for L — 7 in the SF 
phase: from top to bottom, J/U = 0.14, 0.17, and 0.20. 



In the limit L — > oo, one in principle would be able to determine the asymptotic behaviour of the correlation function 
which for large I should behave as 

At the phase boundary, one expects ^oo to diverge resulting in an algebraic decay of the correlation function. 
However, the finite size of the L-mers considered in our calculations precludes extracting useful information about 
this asymptotic behaviour. 

We also observe that in the MI phase the values of Ci{L) for a fixed I depend weakly on L. This is confirmed by 
the lowest curve in Fig. [TUTa) for J/U = 0.1 where we plot all the available data for Ci{L) for L = 3, . . . , 8. This 
shows that the short-range behaviour of the Ci{L) correlation function is not affected significantly by increasing the 
cluster size. The results in the SF phase are shown for L = 7 in Fig. fTUlbV The behaviour shown here is qualitatively 
different from the MI phase. 

In Fig. [TT] we show the correlation length £,{L) obtained by fits of the data to (IMl) as a function of J/U . In the 
MI phase there is a monotonic increase of ^ with J /U , reaching a peak (or cusp) at the phase boundary between the 
MI and SF phases, the position of which of course depends on L. As stated earlier, the variation of ^ with L at a 
fbced valued of J/U is relatively weak in the MI phase; part of this variation is simply a result of fitting data over 
an increasingly larger range of I. However, the L-dependence of ^{L) is much stronger in the SF phase. The distinct 
behaviour of ^(L) between the MI and SF phases persists for all values of n/U. As for the other fitting parameters, we 
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FIG. 11: The correlation length £,{L) vs. J/U for the linear chain of length L = 4, 6 and 8 (bottom to top). The results were 
obtained by fitting the data to (|69|l : the vertical bars indicate the estimated error. 

find that Zq — 1 in the MI phase, reflecting the fact that Ci=o{L) ~ 1; Zq then increases shghtly as Ci=q{L) decreases 
in the SF phase. Finahy, we find that i]{L) is close to 0.4 for most of the data considered. 

Rather similar behaviour is found for the two-dimensional 3x3 cluster in Fig. HJb) in the MI phase. For this 
cluster there are five distinct intersite separations r^j, namely 1, \/2, 2, -s/S, and 2\/2 in units of the lattice constant. 
Due to the finite size of the cluster, there are pairs of sites, such as the nearest-neighbour AB and BC sites, for 
which the correlations are not equivalent. However, the differences are found to be small and we therefore average 
the correlations over all pairs of sites which, by symmetry, would be equivalent in the infinite lattice. These averaged 
values are plotted in Fig. [12] 

The data in Fig. 1121 were fit to a function having the following form 

Cij = e"'''^/« (a -I- b{l - 6ij) cos4%) , (71) 

where ^ is again the correlation length, and the parameters a and b are introduced to capture the anisotropy of the 
lattice correlation function. The variable 9ij is the angle that the vector makes with the a;-axis. The values of Cij 
as determined by (|7T|) are also plotted in Fig. [1^ for the optimal values of the fitting parameters. These values are 
joined by the dotted fines to provide a guide to the eye. It can be seen that the data points fie quite close to these 
lines for all J/U , indicating that the assumed form of the angular dependence does a reasonable job of representing 
the data. The correlation length ^ is plotted in Fig. [T3]as a function of J/U. The increasing trend is similar to that 
found for the linear chain, but the 2D correlation length is significantly smaller than in ID. We emphasize again that 
the correlation length we have extracted represents the short-range behaviour of the correlation function. One would 
need the long-range behaviour in order to determine the critical behaviour of the correlation length at the Ml-SF 
phase boundary. 

Finally, we consider the occupation number distribution for a particular site in the L-mer in the MI phase. The 
probabfiity of finding the configuration in |0), the ground state of /C"' Thus the occupation 

number distribution for, say site 1, is 

P{ni)^ J2 l(W}|0)P- (72) 

712, ••• 

Near the tip of the first Mott lobe of the 3 x 3 L-mer of the 2D square lattice {fi/U = 0.4, J/U = 0.049) we find 
P(0) = 0.022, P(l) = 0.955 and P(2) = 0.023 for the central site of the 3 x 3 i-mer. Even though the average 
occupancy of the sites within the L-mer is exactly n = 1 in the MI phase, the occupancy of a given site does fluctuate 
about its average value; the probability of finding or 2 particles on this site is roughly 2%. 
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FIG. 12: [Colour online] The variation of the averaged pair correlations with separation in units of the lattice constant for the 
two-dimensional 3x3 cluster. The points are the result of the MSMFT calculation for /i/t/ = 0.4 and for J/U — 0.01 (circle), 
0.02 (square), 0.03 (diamond), 0.04 (triangle) and 0.049 (cross). The values of Cij as obtained from the fitting function in (|7ip 
are also plotted as solid points (red) which are joined by dotted lines as a guide to the eye. The MI-SF transition for this value 
of fi/U occurs close to J/U = 0.0492. 
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FIG. 13: The dependence of the correlation length, ^ (in units of the lattice constant), on J/U for the two-dimensional square 
lattice using a 3 x 3 cluster. The results were obtained by fitting the data in Fig. [12] to (|71ll . 
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IV. APPLICATION TO SUPERLATTICES 
A. The Dimer Chain 

We now turn to an analysis of a simple superlattice consisting of a one-dimensional chain in which the on-site 
energies alternate in an ...eaSb^aSb--- fashion. Figure [T^ shows a portion of such a lattice, with a pair of AB sites 
being referred to as a dimer. The BH Hamiltonian for this system can then be expressed as 

jT.i^f~^f)- ^^T.if^f + ^f) + lT.iUAnf{nf 1) + Usnfinf - 1)) 

j 3 j 

- Ji J2 + ^^3) + ' (73) 

3 3 

where the index j labels the jth dimer within the lattice. Here, we have introduced the site operators aj(a|) and 




FIG. 14: Schematic of a general dimer superlattice with alternating site energies ea and eb- The level separation is A. The 
different barrier heights along the chain lead to different intra-dimer (Ji) and inter-dimer ( J2) tunnelling energies. 

bj{bj) for the A and B sites, respectively, in the jth dimer, and the number operators = Ojdj and = b^jbj. 
With EA = —A/2 and eb = A/2, the level separation within a dimer is A. We have also allowed for different 
interaction parameters Ua and Ub for the two kinds of sites and for different intra-dimer (Ji) and inter-dimer (J2) 
hopping (tunnelling) energies. For the optical lattice potential illustrated in Fig. 1141 the asymmetric barrier heights 
lead to Ji > J2. In the following, we consider the optical lattice shown in Fig. 1151 which has inversion symmetry and 
symmetric barrier heights. In this case, Ji = J2 = J. For simplicity, we also assume Ua = Ub = U. This particular 
model was previously studied in Q and p^ . 




FIG. 15: A portion of the dimer superlattice used in the analysis of the MI-SF transition. The inversion symmetry about each 
lattice site results in equal intra- and inter-dimer tunnelling energies J. 

To implement the MSMFT, we partition the lattice into clusters containing Nd dimers, thus defining _L-mers of 
size L — 2Nci. Each L-mer begins with an A site and ends with a B site, and the different L-mers are decoupled 
by introducing the order parameters ipA and tpB- Following the procedure given in IIIBl we obtain the mean- field 
decoupled L-mer Hamiltonian 

IC^'^ =JCl+Vt'^, (74) 

where 




-JY.{% + &]a,) - J E + (75) 
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is the Hamiltonian of an open-ended dimer chain of length L = 2Nd, and 



-J^pB (a\ + ai) - JV-A (^Iv^ + + "iJtpAi^B (76) 



is the mean-field coupling. 

To determine the MI-SF phase boundary we use the energy criterion based on the eigenvalues of the energy matrix 
W. For this example of a two-component order parameter, the expansion of the grand potential has the form given 



in (l58l) with the matrix elements 



^---^^^„f^o({0})-17.({0})' 



and Wba — Wab- The W eigenvalues in this case are given by 



^± - liWAA + Wbb) ± ]J (^l^^^i—^) ' -,- Wl^. (80) 

Since Waa and W^bb are negative definite, we see that w- < 0. On the other hand, since Wab — J for small J, while 
Waa and Wbb are proportional to J^, we see that a;+ remains positive with increasing J up to the point where 

WaaWbb = M^ls, (81) 

when it goes to zero. This equation defines the phase boundary between the MI and the SF phases. We thus see that 
■0 = is a saddle point of the grand potential in the MI phase and that this point turns into a local maximum when 
the SF phase is entered. The stationary point remains a saddle point in the SF phase. 

1. Results for Nd ~ 1 

In this subsection we present results for the case of an L-mer consisting of a single dimer, that is, Nd — 1. This 
case already exhibits all the general features of the phase diagram for a dimer chain. The results we obtain in the 
following subsection show quantitative improvements with increasing but are qualitatively similar. The results we 
find for Nd — 1 are also similar to those found in an earlier investigation [22]. This latter work, however, is a SDMFT 
(even though different order parameters are introduced for the A and B sites) since the dimer Hamiltonian used in 
the perturbative analysis is the sum of independent site Hamiltonians. In contrast, ours is truly a MSMFT, even for 
Nd = 1, since it requires the determination of the eigenstates of the coupled dimer Hamiltonian in ([75)1 . 

The phase diagram for the dimer chain is determined in the three-dimensional parameter space defined by the 
dimensionless variables fl = ii/U , J = J/U and A = A/U. For given values of fi and A, the system is in a MI phase 
for J < Jcr(/i, A). This function defines the phase boundary as a surface in the three-dimensional parameter space. 
As we shall see, this surface consists of sections, each of which terminates on the J = plane and corresponds to 
distinct MI regions. In order to understand the underlying structure of the phase diagram, it is therefore useful to 
first consider the limiting case of J = which defines the base of the Mott lobes. For J = and Nd = 1, we have the 
dimensionless dimer Hamiltonian (we can dispense with the dimer index j in (1751) since there is a single dimer in the 
L-mer) 

}Cl{J = Q)= (-^-li\nA + \ ^-~AnB + \{nA{fiA-l)+nB{nB~l))- (82) 
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FIG. 16: The Af-domains of a single-dimer L-mer at J = 0. The site occupations of the A (lower level) and B (upper level) sites 
is indicated by the number of black dots. As one crosses a downward (upward) sloping line with increasing fJ,/U, the occupation 
of the lower (upper) level increases by one. This figure applies to any number of dimers in the L-mer being considered. 

The occupation number states \nA, tlb) are eigenstates of this Hamiltonian with eigenvalues 

kl^\{nA-'^A^^^ -\{r^A+^^ +\{nB~nB^]^ -\(j^B+^^ , (83) 

where {Ia = fi + A/2 and fiB = A* ^ A/2. The lowest energy is obtained by minimizing with respect to ua and 
UB- Since ua is an integer, the minimizing value of ua is the integer closest to jiA + 1/2, that is, the integer va 
when VA — ^ < fj-A < i^A- This condition defines a region in the /i-A plane where the number of atoms on the A site 
is VA- These regions are shown in Fig. [11] and are bounded by the lines /i = — A/2 for = 0, 1, • • • . A similar 
analysis determines vb , the minimizing value of ub , and the regions where this value applies are bounded by the lines 
/i = + A/2 with vb = 0^ ■ ■ ■ ■ As shown in Fig. I16[ the fl-A plane is thus divided into domains in which the 
number of particles in the dimer is iV = + i^^- These domains are the iV-domains introduced earlier and define the 
base of the distinct MI regions. We observe that each of these regions has specific site occupations. For example, for 
TV = 2, there is a region where the dimer configuration is given by va = vb = ^ and another with va = 2, vb = 0. 
These two configurations are degenerate in energy at a single point in the /i-A plane. They are also degenerate with 
other configurations having a different value of N along the lines bounding the iV-domains. With increasing J, the 
occupation number states are no longer eigenstates of the dimer Hamiltonian but within each Mott lobe the dimer 
configuration is predominantly that indicated in Fig. [12] at the base of the lobe. 

In Fig.[T7]we show a contour plot of the Jcr(/i, A) surface, below which the MI phase exists. The grey scale indicates 
the extent of the MI phase in the J direction which is perpendicular to the plane of the figure. In general, we see that 
the area enclosed by a given contour decreases as J increases, resulting in a dome-like structure of the MI-SF phase 
boundary. The J = contour (black) corresponds to the lines in Fig. (TH] and defines the base of each of the Mott 
lobes. Within each lobe the number of particles in the dimer, {ua) + (fiB), is equal to the integer indicated by the 
configuration in Fig. 1161 However, as we show in more detail later, the average site occupancies (tt-a) and {ub) are 
not constant throughout a given lobe. 

Although Fig. [T7] provides the overall structure of the phase diagram, it is useful to consider various two-dimensional 
cross-sections to visualize the MI and SF regions. A slice through Fig. [T7] at a constant value of J gives one of the 
contours in the figure and reveals the Mott regions as islands in the fl-A plane surrounded by SF regions. Another 
possibility is a vertical slice through Fig. [17] at a fixed value of A. In Fig. [18] we show such a slice at A = 0.75. 
Together with Fig. [T7]it is easy to visualize the variation of the MI regions with variations in A. 

Perhaps a more physical situation is that for a fixed optical lattice, hence fixed values of J and A, but with varying 
U which can be controlled experimentally by means of a Feshbach resonance. A cross-section at fixed J/ A (i.e., J/ A) 
corresponds to a plane containing the fi axis and inclined at some angle with respect to the J and A axes. In Fig. [19] 
we show the intersection of such a plane with the Mott lobes in Fig. [T7] for J/A = 0.05. The superfiuid phase is 
indicated by the 'shaded' region (the vertical black lines) surrounding the various MI phases that occur for different 
values of N. One can see that odd values of N, corresponding to a half-integral average number of particles per 
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FIG. 17: The phase diagram of a dimer chain within the MSMFT using a cluster consisting of a single dimer. The grey scale 
indicates the value of Jcr(/i, A) at which the MI-SF phase boundary occurs. Each region bounded by a black contour (J = 0) 
corresponds to a different Mott lobe with an integer number of atoms N in the dimer. 




J/V 

FIG. 18: A slice through the phase diagram in Fig. [17] at A = 0.75. The average filling per site for each Mott lobe increases 
by 1/2 with increasing /i/(7, starting with n = at the bottom of the figure. The system is a SF to the right of the phase 
boundary. 

site, give rise to lobes which emerge at integral values of /i. These lobes should not be confused with the so-called 
"loopholes" discussed in [Tsl, which arise for a different reason when more complex superlattice structures are 
considered. 

Also shown in Fig. [12] are the various iV-domains which are bounded by the solid curves. These domains evolve 
continuously from those shown in Fig. [Tn]as the J/A = plane is tilted into the J/A = 0.5 plane. As one can see, 
the crossings of the domain boundaries in Fig. [TC] become avoided crossings as a result of the finite value of J. These 
avoided crossings become more pronounced with increasing A since J = 0.05A. 

2. Results for Nd > 1 

We next consider the effect of increasing the number of dimcrs in the L-mer to Nd — 2. The J — phase diagram 
is the same as in Fig. [TC] but the site occupancies shown should be understood to apply to both dimers so that the 
number of particles in each of the A^-domains is actually doubled. The difference from A^^ — 1 becomes apparent for 
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J 7^ and in Fig. [20] we show the phase diagram for J/ A = 0.05. This phase diagram is very similar to the Nd = 1 
phase diagram in Fig. [TO] apart from the appearance of a new Mott region for = 6 corresponding to an average 
fining (n = N/L) of n = 3/2 particles per site. This simply reflects the fact that the Mott regions are stabilized with 
increasing L-mer size as found for the homogeneous systems discussed in mill As increases, Jcr(/i, A) increases 
and the Mott lobes in Fig. [TTl extend to larger values of J. The intersection of the J/ A — 0.05 plane with the Mott 
lobes can thus include additional Mott regions as found in the present example. 
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FIG. 19: The phase diagram of the dirtier chain obtained for an L-mer consisting of a single dimer at J/A = 0.05. The SF 
phase is indicated by the shaded region (the vertical black lines) which surround the various MI regions. The black lines are 
the boundaries of the A'^-domains within which Mott regions with TV particles per dimer are located. As discussed in the text, 
the A''-domains are the regions where the ground state of /C^ contains A'^ particles. 
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FIG. 20: The phase diagram of the dimer chain at J/A = 0.05 for an L-mer containing two dimers. As compared to the single 
dimer phase diagram at in Fig. 1191 we see that the Mott regions are larger, implying the increased stabilization of the Mott 
phase. This is particularly evident in the bottom-right corner of the figure which shows an additional Mott region with filling 
n = 3/2. 

The other qualitative difference in Fig. [20] concerns the iV-domains which arc doubled in number as compared 
to those shown in Fig. [19] Depending on the value of fi the ground state of the Hamiltonian in ([75]) can have 
N = 0, 1, 2, • • • particles and each of these values appears as a distinct A^-domain in Fig. [20] Since L — 2Nd = 4, 
the average number of particles per site in the L-mer is n = N/A. The additional A'^-domains, however, have no 
qualitative effect on the phase diagram and Mott phases only appear at integral and half-integral values of n. We 
have performed some additional calculations for Nd = 3, that is for L-mers with three dimers. As expected, we find 
only slightly expanded MI regions resulting from the increased stabilization of the Mott phases. 
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It is of interest to compare the MSMFT phase diagram with that given in Fig. 17 of Ref. 0] which was obtained 
using the quantum Monte Carlo (QMC) method. These results can effectively be taken as exact. In terms of our 
variables, the data in this figure correspond to ,//A = 0.25, a factor of five larger than the value used in Figs. [19] and 
BUI Before comparing with the data in Ref. [?] , we first note that the extent of the MI regions in the QMC simulations 
for ID systems are enhanced significantly as compared to the MSMFT results (as seen, for example, in Fig. [7]). In 
other words, the Jcr(/i, A) surface of the QMC calculations lies considerably above that of the MSMFT calculations. 
As a result, a direct comparison of the results obtained using the two methods is not particularly meaningful for one 
and the same fixed value of J/A. However, if J/A is scaled in proportion to the relative extent of the MI regions, 
one mi^it expect a correspondence between the two sets of results. This is indeed the case. The results in Fig. 17 
of Ref. IT'I are qualitatively similar to those of Fig. [5T] obtained for J/A = 0.1. One obvious qualitative difference is 
that the MI regions in the MSMFT calculations do not exhibit the pointed shape seen in the QMC results as was 
found for the case of the homogeneous ID lattice (see Fig. [T]). One can infer from this feature that the QMC version 
of Fig. [17] would have cusps running along the top of the Mott lobes. Such cusps are absent in higher dimensions. 
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FIG. 21: The phase diagram of the dimer chain at J/A = 0.1 for an L-mer containing two dimers. The range of fi/U corresponds 
to that in Fig. 17 of Ref. [3] (note that their interaction parameter is U' — U/2). The lowest unshaded region in the figure is a 
n = Mott insulator. 

Finally, it is instructive to consider the variation of the order parameters ipA and ipB and the on-site densities {fiA) 
and {ns) as a function of A along the line fl = 0.5 in Fig. 1201 The order parameters are determined self-consistently 
using the iterative scheme in (j54p . In the MI phase, the order parameters converge to zero and the Hamiltonian ICl 
reduces to the dimer Hamiltonian /C" . Thus in the Mott regions, {fiA) is equal to {nA)oi the value obtained for the 
ground state of K-'j^, with a similar equality for the B site. In the SF regions, (fiA) is no longer equal to {nA)o but the 
latter quantity is still well-defined and it is of interest to compare it with the actual on-site density (n^)- 

In Fig. [21] we plot the on-site densities for the A and B sites terminating the L-mer, as a function of A, for 
J/A = 0.05 and p. = 0.5. We first consider (n^)o and (n_B)o which are shown by the dashed curves in Fig. [H] 
Referring to Fig. [50] we see that the jl = 0.5 line remains in the = 4 domain (n = 1) for A between and 
approximately 3. For this range of A, (n^)o and {nB)o vary continuously with A. For small A, they are both close to 
1 but as A passes through 1, {nA)o increases continuously to about 2 while {nB}o decreases to about 0. This change 
in occupation refiects the change in the dimer state from predominantly |1, 1) to |2, 0) as Fig. [TBI for J/A = would 
suggest. For larger A, the jj, = 0.5 line successively enters iV = 5, 6 and 7 domains. When an A^-domain boundary 
is crossed, the densities {nA)o and (n_B)o calculated for the ground state of JC^ change discontinuously. On entering 
the N = 5 domain, {nA)o jumps approximately to 2.25 (the occupancy of the interior A-site jumps to about 2.75). 
When the A^ = 6 domain is entered, {nA)o jumps to approximately 3. The latter indicates that the K.^ ground state 
now has predominantly the |3,0) configuration. 

The full curves in Fig. [22] show the variation of {ha) and (ns). These quantities generally follow the variation of 
{nA)o and (^5)0, respectively, and, as stated earlier, are equal to the latter in the MI phases where (ipA) and (V's) 
are zero. Unlike {nA)o and (^5)0, however, (n^) and (ns) vary continuously through the SF regions from one MI 
phase to another since in these regions the ground state of JC'^^ is not a number eigenstate. We also find that the 
densities at the interior A and B sites do not deviate much from the densities at the terminating A and B sites of the 
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FIG. 22: A plot of the on-site densities for the terminating A and B sites of an L-mer containing two dimers, as a function 
of A, for the values of J/ A = 0.05 and = 0.5. The densities {ua} and (ns) are obtained from the ground state of the full 
Hamiltonian (|74|) : for comparison, we also show the densities {nA)o and (ns)o corresponding to the ground state of JCl in (|75|) . 
The discontinuities in (71^)0 and {ns)o occur at the boundaries of the TV-domains. On the other hand, (wa) and (ns) vary 
continuously, but show kinks at the boundaries between the MI and SF phases. 
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FIG. 23: A plot of the order parameters, t/ja and ^pB, as a function of A for J/A = 0.05 and jl = 0.5. The order parameters 
are shown for L-mers containing one and two dimers corresponding to the phase diagrams in Figs. \W\ and 1201 

In Fig. [221 we plot the order parameters, tpA and i/'Sj as a function of A for the same values of J/A = 0.05 and 
ft — 0.5. For comparison, we show the order parameters for both the Nd — 1 and Nd ~ 2 L-mers corresponding to 
Figs. [Tni and [201 The /i = 0.5 line passes through two MI regions in the case of Nd = 1 where the order parameters 
are zero, while in the case of Nd = 2, the line passes through three MI regions. The proximity of a MI region in the 
Nd — 1 case leads to the dip in the order parameters seen in the range A ~ 3.5-4. With increasing Nd, the MI phase 
is stabilized and results in the jj, — 0.5 scan passing through a MI region for this range of A. 

In the vicinity of a MI-SF boundary the order parameters are seen to behave as |A — Acr|^^^ which is consistent 
with the mean- field behaviour given by (|32t . We also observe that the value of i/'A is larger than that of V'b, implying 
a larger superfluid fraction on the A-site as compared to the _B-site. This is to be expected given the lower on-site 
energy sa- The condensate density ratio (ipB , however, is somewhat larger than what one might expect on the 
basis of the on-site densities. For example, at A = 1 we have — 0.44 while (?i_B)/(nyi) — 0.33. A similar 

enhancement of the condensate ratio is observed near A = 3 and A = 5. 
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The above results were found for a one-dimensional system, and while one does not expect mean-field theory to 
represent the exact ground state in this case, it is simplest to demonstrate the application of MSMFT to lattices with 
a non-monatomic basis in one dimension. We wish to stress, however, that the qualitative behaviour obtained for the 
dimer chain is also obtained for higher-dimensional lattices. As one example, we have considered the two-dimensional 
honeycomb lattice in which the site energies within a single hexagon alternate ABABAB. The application of MSMFT 
to this system leads to results which are essentially identical to those obtained for the dimer chain. For a value of 
J/A = 0.05 we find a phase diagram similar to that shown in Fig. 1211 the one minor difference is the extent of 
the various Mott lobes. As for the dimer chain problem, the phase diagram for the two-dimensional lattice can be 
understood in terms of the A^-domain structure shown in Fig. 1161 

B. A Superlattice with Loopholes 

In all of the examples we have considered so far, the phase diagrams obtained using SDMFT and MSMFT are 
qualitatively similar. As shown previously by Buonsante et al. p^. [l7j. this is not always the case. In particular, they 
find that certain superlattices exhibit a so-called 'loophole' structure which corresponds to Mott domains that do not 
arise in SDMFT. 

Here we briefly consider an example of this kind. The specific superlattice we study was previously examined using 
SDMFT pJi]. It consists of a four-site superlattice with site energies given by ei = 1.9A, e2 = 0.3A, £3 — 1.3A, 
£4 = O.OA with the same hopping parameter J between nearest-neighbour sites and the same on-site interaction U . 
The phase diagram as determined in SDMFT and MSMFT is shown in Fig. [Ml We observe that Mott domains for 
integral average filling do not arise with SDMFT. However, within MSMFT, these domains are present and have the 
loophole structure referred to above. 
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FIG. 24: Phase diagram for a four-site superlattice which exhibits loopholes. The parameters defining the superlattice are 
given in the text. The shaded region corresponds to the SF phase as determined using MSMFT. The dashed lines show the 
phase boundaries determined using SDMFT. The numbers n labelling each of the Mott domains is the average filling of the 
superlattice. We note that new Mott domains appear at n = 1 and n = 2 within MSMFT. The solid lines show the A^-domains 
for the four-site cluster. Within SDMFT, the A^-domains boundaries are straight lines, parallel to the J/U axis, which start at 
the edges of the Mott domains along the n/lj axis. 

The origin of the integral-n loopholes is associated with degeneracies of different particle configurations in the 
atomic limit (J — > 0). Denoting the /C° ( J = 0) eigenstates as n2, na, 714), the A^ = 4 configurations |0, 2, 0, 2) and 
|0, 1, 1, 2) have the same energy of i?4 = 2.6t7. As the chemical potential is swept through E4 — — 1.3U, the 
particle number in the system jumps by 2. This explains the absence of a n = 1 Mott domain along the J/U — axis. 
As J/U increases, the degeneracy between the two A^ = 4 configurations is lifted within MSMFT and a n = 1 region 
becomes accessible for the formation of a Mott domain. Such a Mott domain does indeed form as seen in Fig. [Ml 
Within SDMFT, the degeneracy persists since /C° (J = 0) is used in the perturbation analysis for all J. As a result, 
the n = 1 Mott domain does not appear. A similar argument applies to the n = 2 Mott domain. Underlying its 
formation is the degeneracy of the |1, 3, 1, 3) and |1, 2, 2, 3) A^ = 8 states in the atomic limit. 

More generally, qualitative differences between SDMFT and MSMFT may be found whenever degeneracies of the 
kind discussed above appear in the ground states of the JCl{J = 0) grand Hamiltonian. Such degeneracies appear in 
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all of the examples of loopholes considered previously [15|, ll7| . They also arise in models dealing with two bosonic 
species 22| and one can expect differences between SDMFT and MSMFT to appear in these cases as well. 



V. CONCLUSIONS 



In this paper we have developed, explained and utilized a mean-field theory of the Mott insulator-superfluid tran- 
sition of bosons moving in optical lattices. In this approach, the lattice is partitioned into clusters which we refer to 
as L-mers. In the Bose-Hubbard model the L-mers are coupled by the hopping Hamiltonian; they can be isolated by 
invoking a mean-field decoupling procedure whereby a superfluid order parameter is introduced for all boundary sites 
of the L-mer. This leads to a mean- field Hamiltonian taking the form shown in Eq. ()12[) . The resulting theory is thus 
a multi-site mean-field theory (MSMFT); it should be contrasted with site-decoupled mean- field theories which allow 
for different order parameters at different sites [l^l but neglect inter-site correlations. 

The ground state of the mean-field Hamiltonian defines a grand potential energy functional of the various order 
parameters. The stationary points of this functional are shown to correspond to self-consistent solutions in that 
the order parameters evaluated at the stationary point coincide with the order parameters defining the mean-field 
Hamiltonian itself. A detailed analysis reveals that the stationary points are, in general, saddle points. As a result, 
they cannot be located by minimizing the energy functional. For weak hopping there is only one self-consistent 
solution, the Mott insulating phase, whereas at larger hopping a second stationary point appears with lower energy, 
corresponding to the superfluid phase. 

The identification of the phase boundaries separating the Mott insulating and superfluid phases can be analyzed 
using perturbation theory. Our work has clarifled the relationship between two different criteria for the determination 
of the phase boundaries, one based on the energy functional itself and the second based on the stability of an iterative 
map. We show that the two criteria are equivalent and can be used interchangeably. 

We have applied our theory to the Bose-Hubbard model for d-dimensional hypercubic lattices, for d — 1 (chains), 
d = 2 (square lattice), and d = 3 (simple cubic lattice). Our results demonstrate the improvements that MSMFT 
affords relative to site-decoupled mean-field theories. This, of course, is what one expects since the theory becomes 
exact in the limit of an infinite L-mer. Specifically, our numerical results make clear that as the size of the L-mer 
and the dimensionality of the system are increased, better agreement between MSMFT and numerical Monte Carlo 
data is obtained. However, mean-field theories typically underestimate the stability of the Mott insulating phase, and 
therefore (for a given chemical potential) underestimate the critical hopping at which this phase becomes unstable 
with respect to the superfluid phase. 

In addition, we have applied our MSMFT to the analysis of one-dimensional superlattices. For the dimer chain, the 
inequivalence of the two sites within the dimer necessitates the introduction of two order parameters. The stationary 
points of the resulting energy functional are found to be saddle points in this case. Apart from the underestimation 
of the critical hopping mentioned above, the phase boundaries we obtain are in qualitative agreement with those 
obtained on the basis of Monte Carlo simulations 7*1 . We also considered one example of a more complex superlattice 
for which the predictions of SDMFT and MSMFT differ qualitatively. Speciflcally, the phase diagram as determined 
using MSMFT exhibits loopholes which are absent in SDMFT. Such qualitative differences are expected to arise when 
degeneracies exist in the JCl{J = 0) ground state energy. 

The MSMFT was used previously in a study of the disordered Bose-Hubbard model ^ISl], and will be applied to 
other, more complicated situations in subsequent work. 



APPENDIX 



A practical limitation of the MSMFT is the size of the Hilbert space needed to diagonalize the mean-fleld Hamiltoni- 
ans /C^ and K.^^^ which grows exponentially with the size of the L-mer. A natural basis is provided by the occupation 
number states ll?!;}) where I = 1, . . . , L and ni = 0, 1, . . . , oo. For a given average number of particles per site, it is 
sufficient to truncate the range of ni at some maximum number rimax which can be varied to ensure convergence of 
the eigenstate calculations. Within this truncated basis, the dimension of the basis is where S — n,nax + 1- 

The perturbative calculations require the ground state of K.^ for N particles and excited states for ± 1 particles. 
The truncated basis discussed above is excessive since it includes states with different total particle numbers. To 
obtain the eigenstates within an A^-particle subspace we therefore retain only those states with the required number 
of particles. We denote the states in this subset as \k), fc = 1, . . . ,Nst, where Ngt is the number of states in the set. 
These states are the occupation number states satisfying X]i=i — ^ with ni < n„iax- 

We now discuss how the size of the basis set can be reduced further with the use of group theory To be specific, 
we consider the 3x3 L-mer illustrated in Fig. [2jb) . The symmetry operations which leave the L-mer invariant define 
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a group Q of order g {g = 8 in this case) with group elements Gm consisting of the identity, three rotations and four 
reflections. The savings provided by group theory stem from the fact that the Hamihonian K-l commutes with the 
group elements Gm and that the ground state |0) of the L-mer belongs to the identity irreducible representation with 
the property Gm,\0) — |0). The reduced energy matrix in ([57)) involves matrix elements of the form X^ieAC*^! 
We observe that 



Gm [aj + aij |0) = Gm J2 + «') G"' G„,|0) 
i£A leA 

^GmJ2{al+ai)G;^'\0) 

leA 

= 5^(a|+a,)|0). (84) 

l£A 

The last step is a consequence of the invariance of the sum of field operators under a symmetry operation. The state 
J2ieA {^l + '^') l'^) t^^^ belongs to the identity representation and as a result, the state \v) must also belong to this 
representation for the matrix element to be finite. It is clear from this discussion that the calculation of the matrix 
in (pS)) cannot be simplified using group theory since the action of Cq, on |0) creates a state that belongs to other 
irreducible representations. 

Since must belong to the identity representation, it can be expanded in terms of states belonging to this 
representation. These states are constructed as follows. 



m — l 



where 1/ ViV^ is an appropriate normalization constant. The number of distinct states formed in this way is Nsym < 
Nst- This reduction in the dimension of the basis set is the advantage afi^orded by the use of group theory. 
Within this symmetrized basis, an eigenstate of K,^ is expanded as 



E ^i''^!^)- (86) 



The expansion coefficients and corresponding eigenvalues r2i,({0}) are then obtained by the diagonalization of the 
Hamiltonian matrix {k\1CI\k!). The reduction of the dimension of the eigenvalue problem from N^t to Nsym is roughly 
a factor 5, the order of the group, and represents a significant computational saving. However, even with the use of 
group theory, the 4x4 L-mer for a 2D square lattice remains a formidable calculation. 
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